# leemos las librerías
library(spatstat)
library(rgdal)
library(sf)
library(ggplot2)

Si ponemos esto, sólo cogeríamos un mapa para poner en la teroía de ppp, pero estoy bloqueada.

Leo Valencia de mapSpain, no me hace falta shapes

# leo Valencia de mamSpain
library(mapSpain)
valencia <- esp_get_munic(munic = "^Valencia$") %>%
  st_transform(25830)

plot(st_geometry(valencia))

# Usamos imagen como mapa de fondo
# tile <- esp_getTiles(valencia, "IDErioja", zoommin = 2)


# Necesitamos un recinto de observación: owin.
val_owin <- as.owin(valencia) # ventana rectangular

DIEGO, DESPUÉS DE MUCHO TOQUETAR, crime_lon TENDRÍA QUE ESTAR MULTIPLICADO POR 10,no?? La cuestión es que estoy atascada. NO se si es por crs o porque soy una inutil :(

## leemos los datos
crime_todos <- read.csv("crime_data_valencia/crime-data-Valencia.csv", header = TRUE, sep = ",", dec = ",")
head(crime_todos)
#>   crime_id crime_date crime_time  crime_type     muni year month week day
#> 1       v1   1/1/2010    4:25:03    Agresion Valencia 2010     1    1   1
#> 2      v10   1/1/2010    0:35:38    Agresion Valencia 2010     1    1   1
#> 3     v100  1/21/2010   19:18:18    Agresion Valencia 2010     1    3  21
#> 4    v1000  4/14/2010   17:12:10    Agresion Valencia 2010     4   15 104
#> 5    v1001  4/15/2010   19:44:02    Agresion Valencia 2010     4   15 105
#> 6    v1002  4/15/2010   21:02:26 Sustraccion Valencia 2010     4   15 105
#>   week_day week_day_name crime_hour    crime_lon   crime_lat    atm_dist
#> 1        5      Thursday          4 -0.382632524 39.46331787 356.2932448
#> 2        5      Thursday          0 -0.390684992 39.43516922 2700.302415
#> 3        4     Wednesday         19 -0.374797136 39.47927475 245.7780012
#> 4        3       Tuesday         17 -0.399208158 39.47982407 891.7241105
#> 5        4     Wednesday         19  -0.35960114  39.4673996 892.4130147
#> 6        4     Wednesday         21 -0.355242282 39.47460938 795.8890846
#>   bank_dist bar_dist cafe_dist industrial_dist market_dist nightclub_dist
#> 1  269.6843 290.3757  240.4669        155.0799    996.8535       458.3678
#> 2   28.2231  38.6437 1221.9344       1010.8148   3270.7254      1387.9167
#> 3  114.4145 232.6933  124.1983        1114.568    438.9616         573.37
#> 4   65.1267 317.7541   253.704       1045.0373    1488.312      1132.8623
#> 5  202.1959   95.727  224.8632        431.1203    576.8191        967.242
#> 6  240.8653 412.9209  115.9555        642.8338    854.4015      1285.3619
#>   police_dist pub_dist restaurant_dist taxi_dist grid_id     grid_lon
#> 1    545.2982 318.2427        150.9095  308.4721     229 -0.384235105
#> 2   2690.2438 829.9902       1373.0596   39.1007     328 -0.389800975
#> 3    901.8223 151.2014        165.2302  445.5541     191 -0.373103365
#> 4   1185.6137 435.2202        258.7675  279.7017     186 -0.400932715
#> 5    746.8792 446.6147         72.0188  458.2727     233 -0.361971625
#> 6    861.3933 308.6051         71.9019 1126.7638     194 -0.356405755
#>      grid_lat
#> 1 39.46532225
#> 2 39.43577975
#> 3 39.47713925
#> 4 39.47713925
#> 5 39.46532225
#> 6 39.47713925
dim(crime_todos)
#> [1] 90247    28

# podemos coger hasta 4 años y tenemos 2010-2020
# cojo, el 2010
mydata_crime <- crime_todos %>% dplyr::filter(year == 2010)


library(sf)
# Objeto sf sin CRS

mydata_crime_sf <- mydata_crime %>%
  st_as_sf(
    coords = c("crime_lon", "crime_lat"),
    crs = st_crs(4326)
  ) %>%
  st_transform(25830)
# Usamos imagen como mapa de fondo
tile <- esp_getTiles(valencia, "IDErioja", zoommin = 2)

ggplot() +
  layer_spatraster(tile) +
  geom_sf(
    data = mydata_crime_sf,
    col = "blue",
    size = 0.3,
    alpha = 0.3
  )

mydata_points <- SpatialPointsDataFrame(
  coords = data.frame(
    x = as.numeric(mydata_crime$crime_lon),
    y = as.numeric(mydata_crime$crime_lat)
  ),
  data = data.frame(rep(0, nrow(mydata_crime)))
)

plot(mydata_points, cex = 0.5, col = "red", pch = "+")


# proj4string(mydata_points) <- CRS("+proj=longlat")
# mydata_points <- st_as_sfc(mydata_points)

Dibujo y veo que el eje x no está bien…

library(maptiles)

valencia_gg <- maptiles::get_tiles(mydata_crime_sf, "Stamen.Watercolor")
plot_tiles(valencia_gg)
plot(st_geometry(mydata_crime_sf), add = TRUE, cex = 0.5, col = "red", pch = "+")

# Creamos el objeto ppp
# Es necesario que todas las coordenadas estén proyectadas en el mismo CRS!
st_crs(mydata_crime_sf) == st_crs(valencia)
#> [1] TRUE

coords <- st_coordinates(mydata_crime_sf)


mydata_ppp <- ppp(
  x = as.numeric(coords[, 1]),
  y = as.numeric(coords[, 2]),
  window = val_owin
)

summary(mydata_ppp)
#> Planar point pattern:  5332 points
#> Average intensity 3.917619e-05 points per square unit
#> 
#> *Pattern contains duplicated points*
#> 
#> Coordinates are given to 1 decimal place
#> i.e. rounded to the nearest multiple of 0.1 units
#> 
#> Window: polygonal boundary
#> 4 separate polygons (no holes)
#>            vertices      area relative.area
#> polygon 1         4   1977930       0.01450
#> polygon 2        82 131953000       0.97000
#> polygon 3         7    958085       0.00704
#> polygon 4         7   1214060       0.00892
#> enclosing rectangle: [720573.4, 735116.6] x [4351274, 4382995] units
#>                      (14540 x 31720 units)
#> Window area = 136103000 square units
#> Fraction of frame area: 0.295
#> 
#> *** 6 illegal points stored in attr(,"rejects") ***

plot(mydata_ppp)

Cálculo de la densidad mediante cuadrantes. Prueba de hipóteis para saber si los puntos se distribuyen al azar o tienen algún patrón.

## Gráfica del objeto de patrones de puntos
# plot(mydata_ppp) no está bien creado
## Resumen del objeto de patrones de puntos.
summary(mydata_ppp)
#> Planar point pattern:  5332 points
#> Average intensity 3.917619e-05 points per square unit
#> 
#> *Pattern contains duplicated points*
#> 
#> Coordinates are given to 1 decimal place
#> i.e. rounded to the nearest multiple of 0.1 units
#> 
#> Window: polygonal boundary
#> 4 separate polygons (no holes)
#>            vertices      area relative.area
#> polygon 1         4   1977930       0.01450
#> polygon 2        82 131953000       0.97000
#> polygon 3         7    958085       0.00704
#> polygon 4         7   1214060       0.00892
#> enclosing rectangle: [720573.4, 735116.6] x [4351274, 4382995] units
#>                      (14540 x 31720 units)
#> Window area = 136103000 square units
#> Fraction of frame area: 0.295
#> 
#> *** 6 illegal points stored in attr(,"rejects") ***

Cálculo de la densidad mediante cuadrantes. Prueba de hipóteis para saber si los puntos se distribuyen al azar o tienen algún patrón.

CC1 <- quadratcount(mydata_ppp,
  nx = 5,
  ny = 5
)
CC1
#> tile
#> Tile row 1, col 1 Tile row 1, col 2 Tile row 1, col 3 Tile row 1, col 4 
#>                 0                12                 0                 0 
#> Tile row 1, col 5 Tile row 2, col 1 Tile row 2, col 2 Tile row 2, col 3 
#>                 0               275              2860              1455 
#> Tile row 2, col 4 Tile row 3, col 1 Tile row 3, col 2 Tile row 3, col 3 
#>               301                 6               278               100 
#> Tile row 3, col 4 Tile row 4, col 2 Tile row 4, col 3 Tile row 4, col 4 
#>                45                 0                 0                 0 
#> Tile row 4, col 5 Tile row 5, col 2 Tile row 5, col 3 Tile row 5, col 4 
#>                 0                 0                 0                 0 
#> Tile row 5, col 5 
#>                 0
plot(mydata_ppp, pch = 1)
plot(CC1, add = TRUE, col = "red")

Estimación de la densidad de patrones de puntos.

plot(density(mydata_ppp), las = 1)
points(mydata_ppp, pch = 19)